The GC skew index: a measure of genomic compositional asymmetry and the degree of replicational selection.

Circular bacterial chromosomes have highly polarized nucleotide composition in the two replichores, and this genomic strand asymmetry can be visualized using GC skew graphs. Here we propose and discuss the GC skew index (GCSI) for the quantification of genomic compositional skew, which combines a normalized measure of fast Fourier transform to capture the shape of the skew graph and Euclidean distance between the two vertices in a cumulative skew graph to represent the degree of skew. We calculated GCSI for all available bacterial genomes, and GCSI correlated well with the visibility of GC skew. This novel index is useful for estimating confidence levels for the prediction of replication origin and terminus by methods based on GC skew and for measuring the strength of replicational selection in a genome.


Introduction
In circular bacterial chromosomes, the replication process starts from a fi nite replication origin (ori) and continues bidirectionally along the two arms (i.e. the replichores) until the replication complex reaches the replication terminus (ter), located directly opposite of ori (Rocha, 2004a;Rocha, 2004b). Replication is obviously the most fundamental and essential process in the cell cycle of bacteria, and replication also exerts genome-wide mutational and selection pressure, shaping genomic polarity with asymmetrically biased nucleotide composition in leading and lagging strands (Lobry and Louarn, 2003;Lobry and Sueoka, 2002). This compositional skew can be easily observed by plotting the normalized excess of guanine (G) over cytosine (C) content in a subgenomic region with sliding windows along the complete genome sequence (Lobry, 1996). Such a GC skew graph segregates the genome into two regions: one with an excess of G over C corresponding to the leading strand, and the other with an excess of C over G corresponding to the lagging strand. Moreover, the shift points of the GC skew graphs are reportedly correlated with the loci of ori and ter (Frank and Lobry, 1999). GC skew is observed in many bacterial species with circular chromosomes, although with varying clarity of the shift points, and GC skew is usually not detectable in symbionts and bacteria with linear chromosomes (Worning et al. 2006) or in archaeal genomes, which employ different machinery for the replication process (Grabowski and Kelman, 2003;Lopez et al. 1999;Myllykallio et al. 2000). GC skew is also observed in local genomic regions primarily introduced by RNA synthesis (Fujimori et al. 2005), but the overall genomic polarity due to replication is present regardless of these local effects, and the GC skew is thus observed in intergenic regions as well as in the third nucleotide positions in codons. Although the underlying causes for GC skew is not completely understood, hydrolytic deamination of cytosine in the leading strand in single-stranded state during replication, is suggested as the major contributing factor (Rocha, 2004b).
Because only a few ori and ter positions had been identifi ed by experimental means, analysis of GC skew was fi rst utilized for the computational prediction of ori and ter positions by examining available genome sequences (Frank and Lobry, 2000). Similar method using nucleotide gradients of T/C and A/G is utilized for the detection of unidirectional replication in mitochondria (Krishnan et al. 2004;Seligmann et al. 2006). To improve the accuracy of prediction, cumulative diagrams are commonly employed to balance out the noise in sequence composition and to eliminate the requirement for window slides (Grigoriev, 1998), coupled with purine and keto excesses and GC skew (Freeman et al. 1998).
However, predictions based on these methods are less accurate in genomes where GC skew cannot be strongly observed (Zawilak et al. 2001). To observe the control of replicational selection on the various genomic properties, genomic compositional skews are also used in conjunction with other genomic features such as the gene orientation (McLean et al. 1998), the distribution of RAG oligomers recognized by the FtsK translocase (Hendrickson and Lawrence, 2006), and the codon bias of genes along the genome (Daubin and Perriere, 2003). To our knowledge, however, no method to quantify the strength of GC skew has been proposed; therefore, it is diffi cult to compare the effects of replicational selection across bacterial genomes.
In this work, we present the GC skew index (GCSI), which quantifi es the strength of GC skew of a given genome by combining Fourier power spectral analysis with the Euclidean distance between the maximum and minimum of the cumulative skew vector. Spectral analysis using fast Fourier transform (FFT) is able to identify the frequency components contributing to a given signal, and it has been applied successfully to the fi eld of bioinformatics (Dodin et al. 2000;Katoh et al. 2002;Yin and Yau, 2005). Because GC skew emerges from the mutational selection in the two replichores, the greatest contributing frequency component of GC skew should be at 1 Hz, with two clear shift points. This observation of a 1-Hz signal combined with the degree of skew calculated by the distance measure between the two vertices of a cumulative skew diagram effectively quantifi es the skew of genomic compositional asymmetry.

Sequences and software
Complete circular chromosomal sequences of 303 bacteria and complete genome sequences of 29 archaeal genomes in GenBank format were selected and obtained from the NCBI RefSeq FTP repository (ftp://ftp.ncbi.nih.gov/genomes/Bacteria/). All analyses were conducted using the G-language Genome Analysis Environment version 1.6.11 (Arakawa et al. 2003;Arakawa and Tomita, 2006). The positional coordinate system for the genomic sequence used in this work was set to originate at 0, unlike that of GenBank, which uses 1 for the position of the fi rst base.

Calculation of GC skew
GC skew was defi ned as the normalized excess of C over G in a given sequence, (C − G)/(C + G), which is calculated with sliding windows along the genome. GC skew is defi ned to be 0 when the amount of C equals that of G. To eliminate the use of window slides, cumulative skew can be calculated as the cumulative sum of the walker graph score at each nucleotide position along the genome, with scores A = 0, T = 0, G = 1, and C = −1. In this work, however, the cumulative GC skew was calculated by taking the cumulative sum of the GC skew in each of the windows, to normalize the cumulative skew strength without it being affected by the length of the genome.

Fast fourier transform
FFT is the computationally optimized derivation of discrete Fourier transform (DFT) for the number of sampling units in the power of two. FFT transforms a given signal in the time domain to reveal the frequency components comprising the input signal. GC skew can be thought of as a signal along the continuous axis of genomic position, which was used in place of the time domain in this work. DFT F(k) of a signal of length N, f (n), where n = 0, 1, …, N − 1, at frequency k was calculated as follows: (1) where i = . The power spectrum PS(k) of F(k) was further defi ned as at each frequency k. In this power spectrum, GC skew shows the greatest contributing component at 1-Hz frequency, corresponding to the two replichores shifting between two regions of opposite polarity as in a sine curve (Arakawa et al. 2007). The Math:FFT module of Perl (http://search.cpan. org/~rkobes/Math-FFT-1.28/FFT.pm) was used for FFT calculation. To level the effects of genome size when comparing the diverse bacterial species, all genome sequences were divided into 4096 windows, and then the GC skew used as the initial signal, the cumulative GC skew, and the power spectra were calculated in these windows. Number of windows must be the power of two for effective FFT calculation, and here 2 12 = 4096 windows were used to take account of the effects of gene positioning, since this window size roughly corresponds the size of genes (about 1kbp) in bacterial genomes. This window size also eliminates other local mutational factors including those within genes, generated by functional requirements in RNA synthesis and translation.

GC skew index
Because cumulative skew should remain around zero under conditions of no strand bias and inversely increase its value in both positive and negative directions where bias is strong, Euclidean distance between the maximal and minimal vertices can be used as a measure of skew. The limitation of this approach and the central challenge for the quantification of genomic compositional skews, however, reside in the mathematical assessment of the skew structure to have exactly two regions physically balanced in length but with opposite polarity of nucleotide content. FFT is a good method for such a purpose, because it is able to reveal the contributing frequency components. Therefore, we used FFT to assess the fi tness of the skew to the replicational selection model and combined this with the Euclidean distance between the two vertices of cumulative skew to calculate the GCSI. The GCSI is defi ned as the normalized average of the Euclidean distance between the two vertices of cumulative skew (dist) and the ratio of spectral strength at 1 Hz and the average strength of spectra in frequency regions 2 Hz or above (SR). Because the replicational selection is the single most dominant factor for GC strand bias, the ratio of spectral strength at 1-Hz frequency and that of all other spectra or their average must be greater than 1. SR was normalized by division with the rounded maximal SR of all bacterial genomes, which we defi ned here as 6000. Likewise, dist was normalized by 600.

Statistical assessment of the signifi cance of GCSI
Signifi cance of the GCSI values is tested using the distribution of GCSI calculated using two sets of randomized data: GCSI calculated using shuffl ed GC skew, where the window order is randomized using the GC skew values calculated with the original genome sequence, and GCSI calculated using shuffl ed genome, where the entire nucleotide sequence of the genome is shuffl ed while conserving the original nucleotide content. Due to calculation costs, statistical test was conducted using 1000 shuffl ed GC skew and 100 shuffl ed genome data sets. Distribution of the resulting GCSI values for the randomized data set was fi rstly tested for its normality using Kolmogorov-Smirnov Lilliefors test, and the signifi cance of the original GCSI value is calculated using the z-score in the distribution of the randomized data set.

Results
To test the applicability of GCSI for the quantifi cation of GC skew strength, we fi rst assessed the correlation between the Euclidean distance of the two vertices of cumulative GC skew, dist, and the Fourier power spectrum ratio, SR, using all genomes ( Fig. 1). The two measures correlated with an R 2 value of 0.6673, showing that the predominance of the 1-Hz frequency component leads to a stronger degree of skew.
Using the measures dist and SR, GCSI was calculated for 304 bacterial genomes; 50 selected species are shown in Table 1 (see supplementary information for comprehensive listings). From the comprehensive list, nine genomes were further selected to illustrate the GC skew graphs plotted with 500 windows at various GCSI values (Fig. 2). As a control, GCSI was also calculated for 29 archaeal genomes, most of which showed no GC skew (Table 2). Because GCSI was normalized by the rounded maximum values of SR and dist, it ranged from 0 to 1. GCSI in bacterial genomes ranged from 0.006 for Gloeobacter violaceus to 0.815 for Clostridium perfringens (mean, 0.207; median, 0.145; SD, 0.173). The majority of archaeal genomes had GCSI <0.05, and the highest GCSI among archaeal genomes (0.122 of Halobacterium sp.) was low compared to those of bacterial genomes. GC skew was not clearly observable in species with GCSI <0.05, but it showed clear shift points when GCSI >0.10. Due to the limited number of iterations, normality test for the statistical assessment using shuffl ed genome sequence did not score well, but that using shuffl ed GC skew passed the test in all genome analyzed. The z-score was generally low and therefore not signifi cant when GCSI <0.05 (especially <0.02), where the GCSI values may not be accurate. On the other hand, GCSI >0.05 scored extremely high z-scores, and therefore these values accurately depict the polarity of the genomes.

Evolutionary Bioinformatics 2007: 3
As can be seen from the GC skew graphs in Figure 2, the degree of skew correlates with GCSI. No skew was observable for G. violaceus and Synechococcus elongatus PCC 7942, with GCSI values of 0.006 and 0.023, respectively, but a gradual rise from negative values to positive values was observed for Synechococcus sp. CC9605, with a GCSI of 0.065, although the skew was not well defi ned. GC skew became visible at a GCSI of 0.098 in Escherichia coli K12, and the clarity was increased in correlation with the GCSI values for scores greater than 1, as represented by the increasing range of the y-axis from ±0.15 at GCSI values around 1 to ±0.4 at a GCSI of 0.815.

Discussion
The nucleotide sequence of a genome is structured and controlled by a myriad of selection pressures, especially in subgenomic regions, as typifi ed by the fact that coding regions are shaped by the essential order and usage of codons. In addition to such requirements in the subgenomic regions, circular bacterial chromosomes experience genome-wide selection through the replication process. The chiral nucleotide composition in the two replication arms is signifi cant; however, with regard to the evolutionary aspects of replicational selection on bacterial chromosomes, no useful method to quantify the degree of genomic compositional asymmetry has been proposed, unlike the wealth of codon bias measures (Suzuki et al. 2005). This lack of indices for genomic compositional skews was likely due to the difficulty of mathematical formulation and detection of the skewing shape of GC skew graphs. To distinguish the degree of skew, we utilized FFT to observe the predominance of the 1-Hz frequency component, which corresponds to the replicational selection on the two replichores, over other frequency components. Combined with the Euclidean distance between the two vertices in cumulative skew graphs, the formulated GCSI captured the strength of GC skew in bacterial chromosomes, as shown by the above results. GCSI scores are diverse even within bacterial genomes with circular chromosomes, ranging from a number Figure 1. Scatter plot of spectral ratio RS against the Euclidean distance between the two vertices in cumulative graph dist. RS measures the goodness-of-fi t of the "shape" of the overall GC skew to be partitioned into two segments corresponding the two replichores, by calculating the relative predominance of the spectral strength of the 1-Hz frequency component over other frequencies upon applying Fast Fourier Transform. dist measures the degree of bias in the leading and lagging strands, by calculating the Euclidean distance between the average GC skew in the two replichores. RS is generally correlated with dist, therefore combination of these two measures as GCSI should correctly represent both the shape of the graph and the degree of skew. of genomes with extremely low values therefore implying the lack of observable GC skew in the genome, to groups of genomes with clear skews as can be seen in Bacilli.
The majority of the archaeal genomes had GCSI <0.05, at which point no noticeable skew is observed even in bacterial genomes. This is also confi rmed by the z-score in the statistical test using randomized data, with low z-scores (therefore implying less signifi cance) when GCSI is less than 0.05. Thus, 0.05 can be employed as a threshold value to determine whether GC skew is present in a genome and therefore whether replicational selection is acting on the organism. Because the GCSI values do not show a Gaussian distribution, however, it should be noted that the indices are not necessarily proportionate with each other. Therefore, GCSI values should not be compared in terms of ratios but in terms of their rank orders. For the direct comparison of quantitative degrees of skew calculated as the ratio of two values, the use of Euclidean distance may be more suitable. However, significant Euclidean distance between the two vertices of cumulative skew may not always result from the polarity exhibited by the GC skew graph; it could also result from local regions of highly biased nucleotide content. Therefore, to ascertain that the skews are controlled by replicational selection, genomes used for such analyses should be position (bp) . GC skew graphs plotted with 500 windows for nine bacteria at different levels of GCSI. GC skew is not observable for the fi rst two species at GCSI <0.05, and becomes evident at GCSI >0.08. At GCSI >0.1, graphs increase their skewness and the shift points and two replichores can be clearly discerned from the graph. Note that the range of Y-axis extends as GCSI values increase. Overall, GCSI correlates with and correctly captures the degree of skew.

Escherichia coli K12
Evolutionary Bioinformatics 2007: 3 Table 1. GCSI, spectral ratio RS, and the Euclidean distance between the two vertices in cumulative graph dist for randomly selected 50 bacterial chromosomes. Signifi cance was calculated using 1) 1000 samples by shuffl ing GC skew windows, and 2) 100 samples by shuffl ing the entire nucleotide sequence of the genome while conserving the nucleotide composition, and the p-value from the normality test and the signifi cance of the original GCSI value using the distribution of randomized samples was given as the z-score. selected beforehand using GCSI or SR at suffi ciently high thresholds (e.g. 0.07 for GCSI and 200 for SR, also noting the z-scores). GCSI would be a useful index for the estimation of confi dence levels for bioinformatics analyses using genomic compositional skews. Predictions of replication origin and terminus by the observation of shift points (i.e. vertices) of cumulative skew diagrams become erroneous when the GC skew is not well defi ned. However, the confi dence level can be easily estimated by taking into account of the magnitude of the GCSI. In this work we have only described the index for GC skew, although the same method is applicable to purine and keto excesses or any other genomic compositional skews, given that the selection is on the two replichores. Similarly, for comparative studies of genomic features related to evolutionary pressures and replication machinery, GCSI can also be used as a measure of replicational selection.